module Mesh2DT0_m
use Geometry_m
use Mesh2D_m, Point => Connect1
use Mesh2DT3_m

type, extends(mesh2D_t) :: Mesh2DT0
    type(Mesh2DT3), pointer :: meshBG
    type(Point), pointer :: elements(:)
contains
    procedure :: initMesh => initMesh2DT0
    procedure :: printNodes => printNodes2DT0
    procedure :: printElements => printElements2DT0
    procedure :: export => export2DT0
    procedure :: exportSolution => exportSolution2DT0
end type

contains
    
subroutine initMesh2DT0(this, mesh2D)
    implicit none
    class(mesh2DT0), intent(out) :: this
    type(mesh2DT3), target, intent(in) :: mesh2D
    
    integer :: nt, nv, i, element(3)
    real*8 :: x(3), y(3), mx, my
    this%meshBG => mesh2D
    nt = mesh2D%nt
    nv = mesh2D%nt
    this%nv = nt; this%nt = nt
    allocate(this%nodes(nv), this%elements(nt))
    
    do i = 1, nt
        element = mesh2D%elements(i)%idx
        x = mesh2D%nodes(element)%xy(1)
        y = mesh2D%nodes(element)%xy(2)
        mx = sum(x) / 3.0D0
        my = sum(y) / 3.0D0
        this%nodes(i) = Node2D((/mx, my/))
        this%elements(i) = Point(i, 0)
    end do
end subroutine

subroutine printNodes2DT0(this)
    implicit none
    class(mesh2DT0), intent(in) :: this
    integer :: i
    do i = 1, this%nv
        write(*, *) this%nodes(i)%xy
    end do
end subroutine

subroutine printElements2DT0(this)
    implicit none
    class(mesh2DT0), intent(in) :: this
    integer :: i
    do i = 1, this%nt
        write(*, *) this%elements(i)
    end do
end subroutine

subroutine export2DT0(this, fileName)
    implicit none
    class(mesh2DT0), intent(in) :: this
    character(len=*), intent(in) :: fileName
    
    call this.meshBG.export(fileName)
end subroutine

subroutine exportSolution2DT0(this, fileName, u1, u2, u3, u4, u5, u6)
    implicit none
    class(Mesh2DT0), intent(in) :: this
    character(len=*), intent(in) :: fileName
	real*8, intent(in) :: u1(this%nv)
    real*8, intent(in), optional :: u2(this%nt), u3(this%nt), u4(this%nt), u5(this%nt), u6(this%nt)

    integer :: element(3), i, j
    type(Mesh2DT3), pointer :: meshBG
    real*8 :: x(3), y(3), x0(3), y0(3), uu(3, 6), t(3, 6)
	open(unit=10, file=fileName)
    meshBG => this%meshBG
	write(10, *) "title = mesh solution file"
    if (present(u6)) then
        write(10, *) 'variables= "x", "y", "u1", "u2", "u3", "u4", "u5", "u6"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(8(F20.10))") x(j), y(j), u1(i), u2(i), u3(i), u4(i), u5(i), u6(i)
            end do
	    end do
    else if (present(u5)) then
        write(10, *) 'variables= "x", "y", "u1", "u2", "u3", "u4", "u5"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(7(F20.10))") x(j), y(j), u1(i), u2(i), u3(i), u4(i), u5(i)
            end do
	    end do
    else if (present(u4)) then
        write(10, *) 'variables= "x", "y", "u1", "u2", "u3", "u4"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(6(F20.10))") x(j), y(j), u1(i), u2(i), u3(i), u4(i)
            end do
	    end do
    else if (present(u3)) then
        write(10, *) 'variables= "x", "y", "u1", "u2", "u3"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(5(F20.10))") x(j), y(j), u1(i), u2(i), u3(i)
            end do
	    end do
    else if (present(u2)) then
        write(10, *) 'variables= "x", "y", "u1", "u2"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(4(F20.10))") x(j), y(j), u1(i), u2(i)
            end do
	    end do
    else
        write(10, *) 'variables= "x", "y", "u1"'
        write(10, *)  "ZONE NODES=", meshBG%nt * 3, ", ELEMENTS=", meshBG%nt, ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE"
        do i = 1, meshBG%nt
            element = meshBG%elements(i)%idx
            x = meshBG%nodes(element)%xy(1)
            y = meshBG%nodes(element)%xy(2)            
            do j = 1, 3
                write(10, "(3(F20.10))") x(j), y(j), u1(i)
            end do
	    end do
    end if

	do i = 1, this%nt
	    write(10, "(3I10)") 3 * i - 2, 3 * i - 1, 3 * i
	end do
	close(10)
end subroutine

end module
